
Job: betamc.std.prg (/home/gov/faculty/ppaolino/parr)
Sat Jun 16 16:15:24 2001, pid=8345, uid=1141(ppaolino)

cheng alg used if cheng=1, johnk otherwise, cheng=       1.0000000 
parm1 =        2.5000000       0.30000000 
parm2 =        1.5000000      -0.30000000 
p =        11.602924 
q =        5.1208106 
true first difference (mean)=      0.24761606 
true first difference (variance)=   -0.0075219259 

descriptive stats for standard beta parameters

       2.6357777       0.29109016 
      0.28527762       0.32201660 
       1.6358210       0.28016751 
     -0.31874352       0.31902720 
 
standard beta mse=     0.010983805 
 
alternative beta mse=     0.010986249 
 
ols mse=     0.011165505 
 
logit trans ols mse=     0.011155427 
 
het ols mse=     0.011328531 
 
logit trans het ols mse=     0.011150665 
 
log trans het ols mse=     0.012347375 
 
first differences for mean and variance for std beta (sim)
      0.24773230      0.041757680 
   -0.0068863529     0.0067891487 
 
first differences for mean and variance for std beta (nosim)
      0.24846459      0.041903750 
   -0.0066786026     0.0065501342 
 
first differences for mean and variance for alt beta (sim)
      0.24638482      0.041560749 
   -0.0066395916     0.0067601514 
 
first differences for mean and variance for alt beta (nosim)
      0.24708399      0.041713976 
   -0.0064282893     0.0065030683 
 
first differences for mean ols
      0.24952259      0.038595977 
 
first differences for mean ols (logit trans)
      0.26105209      0.042181807 
 
first differences for mean and variance for het normal
      0.23913499      0.040209212 
   -0.0066808849     0.0074792507 
 
first differences for mean and variance for logit trans.
      0.25911461      0.044375696 
   -0.0087474123     0.0085211494 
 
first differences for mean ols (log trans)
      0.27668677      0.049617162 
 
efficiency of standard beta vs. alternative beta=       99.534800 
 
efficiency of beta vs. ols=       92.199758 
 
efficiency of beta vs. ols (log trans)=       105.63004 
 
efficiency of standard beta vs. het normal=       98.049346 
 
efficiency of standard beta vs. het normal (log trans)=       109.37750 
 
efficiency of standard beta vs. OLS (log trans)=       137.22354 
 
efficiency of standard beta vs. alternative beta=       99.851977 
 
efficiency of standard beta vs. het normal=       113.96347 
 
efficiency of standard beta vs. het normal (log trans)=       130.35400 
 
standard beta parms' overconfidence=
       105.88138 
       107.64891 
       106.03812 
       108.69907 
 
alternative beta parms' overconfidence=
       104.16123 
       103.13961 
       105.98569 
       105.61067 
 
ols parms' overconfidence=
       90.752193 
       94.967061 
 
ols (log. trans) parms' overconfidence=
       102.87175 
       100.33963 
 
ols (log trans) parms' overconfidence=
       77.883192 
       111.34949 
 
heteroskedastic normal parms' overconfidence=
       101.75726 
       94.721212 
       107.40366 
       101.68602 
 
heteroskedastic log. trans. parms' overconfidence=
       103.53935 
       108.59584 
       111.56752 
       115.38532 
 
standard beta's mean first difference overconfidence
       104.05328 
 
standard beta's variance first difference overconfidence
       96.037618 
 
alternative beta's mean first difference overconfidence
       101.97598 
 
alternative beta's variance first difference overconfidence
       94.075553 
 
check standard beta's percentage within 95% conf int

      0.94500000 
      0.92300000 
      0.93300000 
      0.92200000 
 
check alternative beta's percentage within 95% conf int

      0.92600000 
      0.93900000 
      0.93700000 
      0.92700000 
 
check ols percentage within 95% conf int

      0.96200000 
      0.95500000 
 
check logit trans/ols percentage within 95% conf int

      0.93600000 
      0.94500000 
 
check het ols percentage within 95% conf int

      0.93100000 
      0.95300000 
      0.93500000 
      0.94500000 
 
check logit trans/ols percentage within 95% conf int

      0.93500000 
      0.91400000 
      0.92300000 
      0.89700000 
 
check logit trans/ols percentage within 95% conf int

      0.97700000 
      0.90200000 
 
check percentage of sig parms (standard beta)
       1.0000000 
      0.11700000 
 
check percentage of sig parms (alternate beta)

       1.0000000 
       1.0000000 
       1.0000000 
     0.078000000 
 
check percentage of sig parms (basic ols)

       1.0000000 
       1.0000000 
 
check percentage of sig parms (logit trans ols)

       1.0000000 
       1.0000000 
 
check percentage of sig parms (het ols)

       1.0000000 
       1.0000000 
       1.0000000 
      0.16600000 
 
check percentage of sig parms (het ols - logit trans)

       1.0000000 
       1.0000000 
      0.99200000 
      0.13900000 
 
check percentage of sig parms (ols -- log tranform)

       1.0000000 
       1.0000000 
 
descriptive stats for y's for successful runs

      0.68098971 
      0.16495778 
      0.30153999 
      0.93441118 
descriptive stats for y's for failed runs
       0.0000000 
descriptive stats for x's for successful runs

     -0.30006816 
      0.98784384 
      -2.1863198 
       1.6917737 
descriptive stats for x's for failed runs
       0.0000000 
number of failed runs
       0.0000000 

code used to generate the above results:

library cml;
#include cml.ext;
CMLset;

fails=0;
crash=0;
ystatf=0;
xstatf=0;

/* load x matrix based upon number of observations */

@ xd = loadd("/home/gov/faculty/ppaolino/betamc/xdata_n500"); @
x = ones(30,1)~rndn(30,1);
n = rows(x);

/* produce matrices for simulated first differences */

meanx=meanc(x); @ create vector of means for the indep vars @
stdx=stdc(x);   @ create vector of stdev for the indep vars @
minx=minc(x);
maxx=maxc(x);

/* start obtaining the first difference vectors (means) */

sw=zeros(1,cols(x)-1)|eye(cols(x)-1);
xl=meanx-stdx.*sw;           @ these are kXk-1 vectors @
xh=meanx+stdx.*sw;       
xmin=meanx.*(1-sw)+minx.*sw;  
xmax=meanx.*(1-sw)+maxx.*sw;

/* start creating Monte Carlo data */
k = cols(x);
parm1 = {2.5, 0.3};
parm2 = {1.5, -0.3};
cheng=1;
print "cheng alg used if cheng=1, johnk otherwise, cheng=";;cheng; 
p = exp(x*parm1);
q = exp(x*parm2);

/* determine true first differences */

/* get  alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alt=exp(xl'*parm1);                          @ nxk-1 matrix @
         blt=exp(xl'*parm2);
         aht=exp(xh'*parm1);
         bht=exp(xh'*parm2);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlt=alt./(alt+blt);                    @ nxk-1 matrix @
         meanht=aht./(aht+bht);
         varlt=(alt.*blt)./(((alt+blt)^2).*(alt+blt+1));
         varht=(aht.*bht)./(((aht+bht)^2).*(aht+bht+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmt=meanht-meanlt;
         fdvt=varht-varlt;

/* print data parameters for simulation */

print "parm1 = ";; parm1';
print "parm2 = ";; parm2';
print "p = ";; meanc(exp(x*parm1));
print "q = ";; meanc(exp(x*parm2));
print "true first difference (mean)=";; fdmt;
print "true first difference (variance)=";; fdvt;

/* start Monte Carlo trials */
trial=1;
do while trial<=1000; 

print "MC trial=";; trial;

   if trial==1;
      ret=0;
      if fails==0;
         bt=0; 
      endif;
   endif;

   if fails==0;
   success=0;

/* comment out random number generator if p,q<1 */
if cheng==1;
/* create J&K (Chen) alpha */
ralpha=p+q;

/* create J&K (Chen) beta */
p1=p.>1;
q1=q.>1;
rbeta=real((p1.*q1).*(sqrt((ralpha-2)./(2.*p.*q-ralpha))) +
      (1-p1.*q1).*(1/(minc((p'|q')))));

/* create J&K (Chen) gamma */
rgamma=p+1/rbeta;

/* step 1 in J&K (Chen) RNG */
u1=rndu(n,1);
u2=rndu(n,1);
v=rbeta.*ln(u1./(1-u1));
w=p.*exp(v);

/* step 2 in J&K (Chen) RNG */

done=0;
      do while done<n;
         d=(ralpha.*ln(ralpha./(q+w))+rgamma.*v-1.3862944).<ln((u1^2).*u2); 
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v=(1-d).*v + d.*rbeta.*ln(u1./(1-u1));
         w=(1-d).*w + d.*p.*exp(v);
         done=n-sumc(d);
      endo;

/* step 3 in J&K (Chen) RNG */
y=w./(q+w); 
endif;

/* discretize y 

y1=yc.>0;
y2=yc.>0.14;
y3=yc.>0.28;
y4=yc.>0.42;
y5=yc.>0.56;
y6=yc.>0.70;
y7=yc.>0.84;

y=(y1+y2+y3+y4+y5+y6+y7)/8; */

/* comment out random number generator if p or q>1 */
if cheng/=1;
/* step 1 in Johnk's RNG */

u1=rndu(n,1);
u2=rndu(n,1);

v1=u1^(1/p);
v2=u2^(1/q);

/* step 2 in Johnk's RNG */

w=v1+v2;

/* step 3 in Johnk's RNG */

done=0;
      do while done<n;
         d=w.>1;
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v1=(1-d).*v1 + d.*u1^(1/p);
         v2=(1-d).*v2 + d.*u2^(1/q);
         w=v1+v2;
         done=n-sumc(d);
      endo;

/* step 4 in Johnk's RNG */
y=v1./w;
endif;

/* cheat function to re-set Ys that equal 1 --- would choke llik function
   otherwise */
d=y.==1;
print "cheat=";;sumc(d);
y=(1-d).*y + d.*.99999;  


/* return the monte carlo data set */
    data = y~x[.,2:cols(x)];

   endif;

   if fails==0;

      do while success==0 and fails<=30;

/* maximum likelihood function for estimating parameters alpha and beta */
__title=" ** A Beta Model, Analytic Solution (Standard Estimation) **";

stval=1.0|ones(k-1,1)*0.5|1.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

       if fails>=1 and fails<=40;
          stval=(bt')/2;
          print "retrying run";
       endif;

print "mean y, max y, min y, std y";
meanc(y)~maxc(y)~minc(y)~stdc(y);

proc psi(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;

proc gradproc(theta,ds);
    local y,x,z,be,h,xb,p,q,per1,gr1,gr2,grad;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    per1=p+q;
    gr1= psi(p+q).*x.*p-psi(p).*x.*p+x.*p.*ln(y);
    gr2= psi(p+q).*x.*q-psi(q).*x.*q+x.*q.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);
    endp;

proc loglik(theta,ds);
    local y,x,z,be,h,xb,p,q,l1,l2,l3,l4,l5,llik;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradproc;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400;

   if fails.>=5 and fails.<=10;
     _cml_DirTol=.000001;
   endif;
   if fails.>10 and fails.<=15;
     _cml_DirTol=.0000001;
   endif;
   if fails.>15 and fails.<=20;
     _cml_DirTol=.00000001;
   endif;
   if fails.>20;
     _cml_DirTol=.00001;
   endif;

    {b,logl,g,vc,ret}=CMLprt(CML(data,0,&loglik,stval));

/* check that parameters are reasonable */
   bt=(b');

/* determine success or failure of the trial */
   cor=vc[1,1]/vc[1,1];
         if ret>=1 or cor/=1;
            fails=fails+1;
         elseif ret==0 and cor==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

   if fails==0;
      success=0; bta=0;
      do while success==0 and fails<=30;

CMLset;
/* Alternative beta estimation */
__title=" ** A Beta Model, Numerical Solution (Alternative Estimation) **";

stvala=0.0|ones(k-1,1)|2.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvala=(bta')/2;
         print "retrying run";
      endif;

proc gradpa(theta,ds);
    local y,one,x,z,be,h,xb,zh,zh1,e,per1,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    zh1=zh+1;
    e=xb./(1+xb);
    per1=(-x.*e+x.*(e^2)).*zh1;

    gr1=psi(e.*zh1+(1-e).*zh1-1).*(x.*e.*zh1-(e^2).*zh1.*x+per1)-
        psi(e.*zh1-e).*(x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x)-
        psi((1-e).*zh1-1+e).*(per1+x.*e-(e^2).*x)+
        (x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x).*ln(y)+
        (per1+x.*e-(e^2).*x).*ln(1-y);
    gr2=psi(e.*zh1+(1-e).*zh1-1).*(e.*z.*zh+(1-e).*z.*zh)-
        (psi(e.*zh1-e).*xb.*z.*zh)./(1+xb)-
        psi((1-e).*zh1-1+e).*(1-e).*z.*zh+xb.*z.*zh.*ln(y)./(1+xb)+
        (1-e).*z.*zh.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglika(theta,ds);
    local y,one,x,z,be,h,xb,zh,e,v,p,q,llik,ds1;
		     
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    e=xb./(1+xb);
    v=e.*(1-e).*(1/(zh+1));
    p=((e^2).*(1-e))./v-e;
    q=(e.*(1-e)^2)./v-(1-e);
    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradpa;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400; 
    {ba,logl,g,vca,ret}=CMLprt(CML(data,0,&loglika,stvala));

/* check that parameters are reasonable */
   bta=(ba');

/* determine success or failure of the trial */
   cora=vca[1,1]/vca[1,1];
         if ret>=1 or cora/=1;
            fails=fails+1;
         elseif ret==0 and cora==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

/* linear heteroskedastic estimation */

   if fails==0;
      success=0; bth=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvalh=inv(x'*x)*x'*y|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=10;
         stvalh=(bth')/2;
         print "retrying run";
      endif;  

proc gradph(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikh(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ds[.,1];    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradph;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {blh,logl,g,vch,ret}=CMLprt(CML(data,0,&loglikh,stvalh));

/* check that Monte carlo parameters are reasonable */
bth=(blh');

/* determine success or failure of the trial */
   corh=vch[1,1]/vch[1,1];
       if ret>=1 or corh/=1;
          fails=fails+1;
          elseif ret==0 and corh==1;
          success=1; fails=0;
       endif;
    endo;
endif;

/* linear heteroskedastic estimation (with logit tranformation) */

   if fails==0;
      success=0; btl=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvall=inv(x'*x)*x'*(ln(y./(1-y)))|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvall=(btl')/2;
         print "retrying run";
      endif;  

proc gradpl(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ln(ds[.,1]./(1-ds[.,1]));
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikl(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ln(ds[.,1]./(1-ds[.,1]));    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradpl;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {bll,logl,g,vcl,ret}=CMLprt(CML(data,0,&loglikl,stvall));

/* check that Monte carlo parameters are reasonable */
btl=(bll');

/* determine success or failure of the trial */
   corl=vcl[1,1]/vcl[1,1];
      if ret>=1 or corl/=1;
         fails=fails+1;
      elseif ret==0 and corl==1;
         success=1; fails=0;
      endif;
   endo;
endif;

/* if trial has failed 30 times */

   if fails==30;
      print "run failed after 30 attempts";
      crash=crash+1;

      if ystatf==0;
         ystatf=(meanc(y)~stdc(y)~minc(y)~maxc(y));
      elseif ystatf/=0;
         ystatf=ystatf|(meanc(y)~stdc(y)~minc(y)~maxc(y));
      endif;

      if xstatf==0;
         xstatf=(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|maxc(x[.,2]));
      elseif xstatf/=0;
         xstatf=xstatf~(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|
                 maxc(x[.,2]));
      endif;
      fails=0;
      clear y;
   endif;

/* collect data for a successful trial */
   if success==1;

/* collect y stats for a successful trial */
      if trial==1;
         ystats=meanc(y)~stdc(y)~minc(y)~maxc(y);
         xstats=(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=b;
         sesum=sqrt(diag(vc));
         bsuma=ba;
         sesuma=sqrt(diag(vca));
         bsumh=blh;
         sesumh=sqrt(diag(vch));
         bsuml=bll;
         sesuml=sqrt(diag(vcl));
      endif;
      if trial>=2;
         ystats=ystats|(meanc(y)~stdc(y)~minc(y)~maxc(y));
         xstats=xstats~(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=bsum~b;
         sesum=sesum~(sqrt(diag(vc)));
         bsuma=bsuma~ba;
         sesuma=sesuma~sqrt(diag(vca));
         bsumh=bsumh~blh;
         sesumh=sesumh~sqrt(diag(vch));
         bsuml=bsuml~bll;
         sesuml=sesuml~sqrt(diag(vcl));
      endif;

/* vectors of first difference values for the j^th independent var */
         bns=b[1:cols(x),.];                @ nxk matrix @
         hns=b[cols(x)+1:rows(b),.];   @ nxk matrix @

/* calculate model mse */

         ahat=exp(x*bns);
	 bhat=exp(x*hns);
         yhatns=ahat./(ahat+bhat);
	 bnsmse=sumc((y-yhatns)^2)/n;
	 
/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alns=exp(xl'*bns);                          @ nxk-1 matrix @
         blns=exp(xl'*hns);
         ahns=exp(xh'*bns);
         bhns=exp(xh'*hns);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlns=alns./(alns+blns);                    @ nxk-1 matrix @
         meanhns=ahns./(ahns+bhns);
         varlns=(alns.*blns)./(((alns+blns)^2).*(alns+blns+1));
         varhns=(ahns.*bhns)./(((ahns+bhns)^2).*(ahns+bhns+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmns=meanhns-meanlns;
         fdvns=varhns-varlns;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmns=fdmns;  @ k-1X1 matrix for mean @
            totvns=fdvns;  @ k-1X1 matrix for variance @
	    totbmse=bnsmse;
         endif;

         if trial>=2 and ret==0; 
            totmns=totmns~fdmns;  @ k-1X(trial) matrix for mean @
            totvns=totvns~fdvns;  @ k-1X(trial) matrix for variance @
	    totbmse=totbmse|bnsmse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasim=(rndmn(b,vc,1000));  @ nxk matrix, where n=# of draws @
         be=betasim[.,1:cols(x)];                @ nxk matrix @
         h=betasim[.,cols(x)+1:cols(betasim)];   @ nxk matrix @

/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         al=exp(be*xl);                          @ nxk-1 matrix @
         bl=exp(h*xl);
         ah=exp(be*xh);
         bh=exp(h*xh);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlow=al./(al+bl);                    @ nxk-1 matrix @
         meanhigh=ah./(ah+bh);
         varlow=(al.*bl)./(((al+bl)^2).*(al+bl+1));
         varhigh=(ah.*bh)./(((ah+bh)^2).*(ah+bh+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdm=meanhigh-meanlow;
         fdv=varhigh-varlow;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totm=meanc(fdm);  @ k-1X1 matrix for mean @
            totv=meanc(fdv);  @ k-1X1 matrix for variance @
            totms=stdc(fdm);  @ k-1X1 matrix for mean @
            totvs=stdc(fdv);  @ k-1X1 matrix for variance @
         endif;

         if trial>=2 and ret==0; 
            totm=totm~meanc(fdm);  @ k-1X(trial) matrix for mean @
            totv=totv~meanc(fdv);  @ k-1X(trial) matrix for variance @
            totms=totms~stdc(fdm);  @ k-1X(trial) matrix for mean @
            totvs=totvs~stdc(fdv);  @ k-1X(trial) matrix for variance @
         endif;

/* get relevant stats for alternative method of beta estimation */
/* vectors of first difference values for the j^th independent var */
         bnsa=ba[1:cols(x),.];                @ nxk matrix @
         hnsa=ba[cols(x)+1:rows(ba),.];   @ nxk matrix @

/* calculate mse for alternative beta estimation */
  
         yhatnsa=exp(x*bnsa)./(1+exp(x*bnsa));
 	 bnsamse=sumc((y-yhatnsa)^2)/n;
	 
/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlnsa=(exp(xl'*bnsa)./(1+exp(xl'*bnsa)))';   @ nxk-1 matrix @
         meanhnsa=(exp(xh'*bnsa)./(1+exp(xh'*bnsa)))';
         dlsa=(1/(exp(xl'*hnsa)+1))';
         dhsa=(1/(exp(xh'*hnsa)+1))';
         varlnsa=meanlnsa.*(1-meanlnsa).*dlsa;
         varhnsa=meanhnsa.*(1-meanhnsa).*dhsa;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmnsa=meanhnsa-meanlnsa;
         fdvnsa=varhnsa-varlnsa;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmnsa=fdmnsa;  @ k-1X1 matrix for mean @
            totvnsa=fdvnsa;  @ k-1X1 matrix for variance @
	    totbamse=bnsamse;
         endif;

         if trial>=2 and ret==0; 
            totmnsa=totmnsa~fdmnsa;  @ k-1X(trial) matrix for mean @
            totvnsa=totvnsa~fdvnsa;  @ k-1X(trial) matrix for var @
	    totbamse=totbamse|bnsamse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasima=(rndmn(ba,vca,1000));  @ nxk matrix, where n=# of draws @
         bea=betasima[.,1:cols(x)];                @ nxk matrix @
         ha=betasima[.,cols(x)+1:cols(betasima)];   @ nxk matrix @

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanla=(exp(bea*xl)./(1+exp(bea*xl)));   @ nxk-1 matrix @
         meanha=(exp(bea*xh)./(1+exp(bea*xh)));
         dl=(1/(exp(ha*xl)+1));
         dh=(1/(exp(ha*xh)+1));
         varla=meanla.*(1-meanla).*dl;
         varha=meanha.*(1-meanha).*dh;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdma=meanha-meanla;
         fdva=varha-varla;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totma=meanc(fdma);  @ k-1X1 matrix for mean @
            totva=meanc(fdva);  @ k-1X1 matrix for variance @
            totmsa=stdc(fdma);
            totvsa=stdc(fdva);
         endif;

         if trial>=2 and ret==0; 
            totma=totma~meanc(fdma);  @ k-1X(trial) matrix for mean @
            totva=totva~meanc(fdva);  @ k-1X(trial) matrix for variance @
            totmsa=totmsa~stdc(fdma);
            totvsa=totvsa~stdc(fdva);
         endif;

/* set parameter vectors (linear heteroskeadastic) */
         bhh=blh[1:cols(x),1];
         hh=trimr(blh,rows(bhh),0);

/* calculate mse */

         yhatlh=x*bhh;
	 lhmse=sumc((y-yhatlh)^2)/n;
	 
/* get predicted means */
         elh=xl'*bhh;
         ehh=xh'*bhh;

/* get predicted variances */
         vlh=exp(xl'*hh);
         vhh=exp(xh'*hh);

/* produce first differences */
         if trial==1 and ret==0;
         totfdmh=ehh-elh;
         totfdvh=vhh-vlh;
	 totlhmse=lhmse;
         elseif trial>=2 and ret==0;
         totfdmh=totfdmh~(ehh-elh);
         totfdvh=totfdvh~(vhh-vlh);
	 totlhmse=totlhmse|lhmse;
         endif;

/* set parameter vectors (linear logit tranformation) */
         blo=bll[1:cols(x),1];
         hl=trimr(bll,rows(blo),0);
	 
/* calculate mse for linear het log trans */
  
         yhatlhlt=exp(x*blo)./(1+exp(x*blo));
	 lhltmse=sumc((y-yhatlhlt)^2)/n;

/* get predicted means */
         eol=exp(xl'*blo)./(1+exp(xl'*blo));
         eoh=exp(xh'*blo)./(1+exp(xh'*blo));

/* get predicted variances */
         vol=((exp(xl'*blo)./(1+exp(xl'*blo))^2)^2).*exp(xl'*hl);
         voh=((exp(xh'*blo)./(1+exp(xh'*blo))^2)^2).*exp(xh'*hl);

/* produce first differences */
         fdml=eoh-eol;
         fdvl=voh-vol;

/* produce first differences */
         if trial==1 and ret==0;
         totfdml=fdml;
         totfdvl=fdvl;
	 totthmse=lhltmse;
         elseif trial>=2 and ret==0;
         totfdml=totfdml~fdml;
         totfdvl=totfdvl~fdvl;
	 totthmse=totthmse|lhltmse;
         endif;

/* run ols */

{ vnam,m,bols,stb,vcols,stderr,sigma,cx,rsq,resid,dbw } = ols(0,y,x);

/* calculate mse */

      yhat=x*bols;
      olsmse=sumc((y-yhat)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsum=bols;
         selsum=stderr;
      endif;
      if trial>=2;
         blsum=blsum~bols;
         selsum=selsum~stderr;
      endif;

/* obtain the first differences */

olsl=xl'*bols;
olsh=xh'*bols;
fdols=olsh-olsl;

      if trial==1;
         fdolsum=fdols;
	 totomse=olsmse;
      endif;
      if trial>=2;
         fdolsum=fdolsum~fdols;
	 totomse=totomse|olsmse;
      endif;

/* run ols (with logit trans) */

yl=ln(y./(1-y));

{ vnam,m,bolsl,stb,vcols,stderrl,sigma,cx,rsq,resid,dbw } = ols(0,yl,x);

/* calculate mse */

      yhatlt=exp(x*bolsl)./(1+exp(x*bolsl));
      ltmse=sumc((y-yhatlt)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumo=bolsl;
         selsumo=stderrl;
      endif;
      if trial>=2;
         blsumo=blsumo~bolsl;
         selsumo=selsumo~stderrl;
      endif;

/* obtain the first differences */

olslo=exp(xl'*bolsl)./(1+exp(xl'*bolsl));
olsho=exp(xh'*bolsl)./(1+exp(xh'*bolsl));
fdolsl=olsho-olslo;

      if trial==1;
         fdolsuml=fdolsl;
	 tototmse=ltmse;
      endif;
      if trial>=2;
         fdolsuml=fdolsuml~fdolsl;
	 tototmse=tototmse|ltmse;
      endif;

/* run ols (with log trans) */

ly=ln(y+.01);

{ vnam,m,bolsll,stb,vcols,stderrll,sigma,cx,rsq,resid,dbw } = ols(0,ly,x);

/* calculate mse */

      yhatll=exp(x*bolsll)-0.01;
      ltmsel=sumc((y-yhatll)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumol=bolsll;
         selsumol=stderrll;
      endif;
      if trial>=2;
         blsumol=blsumol~bolsll;
         selsumol=selsumol~stderrll;
      endif;

/* obtain the first differences */

olslol=exp(xl'*bolsll)-0.01;
olshol=exp(xh'*bolsll)-0.01;
fdolsll=olshol-olslol;

      if trial==1;
         fdsumll=fdolsll;
	 totmsell=ltmsel;
      endif;
      if trial>=2;
         fdsumll=fdsumll~fdolsll;
	 totmsell=totmsell|ltmsel;
      endif;

/* reset fails and add a trial for a successful trial */
      fails=0;
      trial=trial+1; 
      clear y;

   endif;

endo;

/* check that parameters are accurate */
print "descriptive stats for standard beta parameters";
meanc(bsum')~stdc(bsum'); 
print " ";
/* check the mse for the std beta */
  
print "standard beta mse=";;meanc(totbmse);
print " ";
print "alternative beta mse=";;meanc(totbamse);
print " ";
print "ols mse=";;meanc(totomse);
print " ";
print "logit trans ols mse=";;meanc(tototmse);
print " ";
print "het ols mse=";;meanc(totlhmse);
print " ";
print "logit trans het ols mse=";;meanc(totthmse);
print " ";
print "log trans het ols mse=";;meanc(totmsell);
print " ";

/* check mean and std. dev of first differences */

print"first differences for mean and variance for std beta (sim)";
meanc(totm')~stdc(totm');
meanc(totv')~stdc(totv');
print " ";

print "first differences for mean and variance for std beta (nosim)";
meanc(totmns')~stdc(totmns');
meanc(totvns')~stdc(totvns');
print " ";

print "first differences for mean and variance for alt beta (sim)";
meanc(totma')~stdc(totma');
meanc(totva')~stdc(totva');
print " ";

print "first differences for mean and variance for alt beta (nosim)";
meanc(totmnsa')~stdc(totmnsa');
meanc(totvnsa')~stdc(totvnsa');
print " ";

/* check for mean and std. dev. of first difference for OLS */
print"first differences for mean ols";
meanc(fdolsum')~stdc(fdolsum');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (logit trans)";
meanc(fdolsuml')~stdc(fdolsuml');
print " ";

print "first differences for mean and variance for het normal";
meanc(totfdmh')~stdc(totfdmh');
meanc(totfdvh')~stdc(totfdvh');
print " ";

print "first differences for mean and variance for logit trans.";
meanc(totfdml')~stdc(totfdml');
meanc(totfdvl')~stdc(totfdvl');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (log trans)";
meanc(fdsumll')~stdc(fdsumll');
print " ";

/* efficiency of expected values */
/* compare efficiency of standard beta with alternative beta */

effnba=((totmnsa-fdmt)^2);
effdba=((totmns-fdmt)^2);
effalt=100*(sqrt(sumc(effnba')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. alternative beta=";;effalt;
print " ";

/* compare efficiency of standard beta with ols */

effno=((fdolsum-fdmt)^2);
effdb=((totmns-fdmt)^2);
effols=100*(sqrt(sumc(effno')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols=";;effols;
print " ";

/* compare efficiency of standard beta with ols */

effnl=((fdolsuml-fdmt)^2);
effdb=((totmns-fdmt)^2);
effolsl=100*(sqrt(sumc(effnl')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols (log trans)=";;effolsl;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbh=((totfdmh-fdmt)^2);
effdba=((totmns-fdmt)^2);
effhet=100*(sqrt(sumc(effnbh')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. het normal=";;effhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbll=((totfdml-fdmt)^2);
effdbal=((totmns-fdmt)^2);
efflog=100*(sqrt(sumc(effnbll')))./(sqrt(sumc(effdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;efflog;
print " ";

/* compare efficiency of standard beta with log trans normal */

effnbl=((fdsumll-fdmt)^2);
effdba=((totmns-fdmt)^2);
efflogl=100*(sqrt(sumc(effnbl')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. OLS (log trans)=";;efflogl;
print " ";

/* efficiency for variance term */
/* compare efficiency of standard beta with alternative beta */

veffnba=((totvnsa-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffalt=100*(sqrt(sumc(veffnba')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. alternative beta=";;veffalt;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbh=((totfdvh-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffhet=100*(sqrt(sumc(veffnbh')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. het normal=";;veffhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbll=((totfdvl-fdvt)^2);
veffdbal=((totvns-fdvt)^2);
vefflog=100*(sqrt(sumc(veffnbll')))./(sqrt(sumc(veffdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;vefflog;
print " ";

/* check standard beta parameters' overconfidence */

npoca=(bsum-meanc(bsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesum^2;
dpoc=sqrt(sumc(dpoca'));
ocp=100*(npoc./dpoc);
print "standard beta parms' overconfidence=";;ocp;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check alternative beta parameters' overconfidence */

npoca=(bsuma-meanc(bsuma'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuma^2;
dpoc=sqrt(sumc(dpoca'));
ocpa=100*(npoc./dpoc);
print "alternative beta parms' overconfidence=";;ocpa;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols parameters' overconfidence */

npoca=(blsum-meanc(blsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsum^2;
dpoc=sqrt(sumc(dpoca'));
oco=100*(npoc./dpoc);
print "ols parms' overconfidence=";;oco;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log. trans) parameters' overconfidence */

npoca=(blsumo-meanc(blsumo'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumo^2;
dpoc=sqrt(sumc(dpoca'));
ocol=100*(npoc./dpoc);
print "ols (log. trans) parms' overconfidence=";;ocol;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log trans) parameters' overconfidence */

npoca=(blsumol-meanc(blsumol'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumol^2;
dpoc=sqrt(sumc(dpoca'));
ocoll=100*(npoc./dpoc);
print "ols (log trans) parms' overconfidence=";;ocoll;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic normal parameters' overconfidence */

npoca=(bsumh-meanc(bsumh'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesumh^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic normal parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic logit transform parameters' overconfidence */

npoca=(bsuml-meanc(bsuml'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuml^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic log. trans. parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check overconfidence for the first differences (mean) */

print "standard beta's mean first difference overconfidence";
nmoca=(totm-meanc(totm'))^2;
nmoc=sqrt(sumc(nmoca'));
dmoca=totms^2;
dmoc=sqrt(sumc(dmoca'));
ocm=100*(nmoc./dmoc);
ocm;
print " ";

/* check overconfidence for the first differences (variance) */

print "standard beta's variance first difference overconfidence";
nvoca=(totv-meanc(totv'))^2;
nvoc=sqrt(sumc(nvoca'));
dvoca=totvs^2;
dvoc=sqrt(sumc(dvoca'));
ocv=100*(nvoc./dvoc);
ocv;
print " ";

/* check overconfidence for the first differences (mean) */

print "alternative beta's mean first difference overconfidence";
nmocaa=(totma-meanc(totma'))^2;
nmoca=sqrt(sumc(nmocaa'));
dmocaa=totmsa^2;
dmoca=sqrt(sumc(dmocaa'));
ocma=100*(nmoca./dmoca);
ocma;
print " ";

/* check overconfidence for the first differences (variance) */

print "alternative beta's variance first difference overconfidence";
nvocaa=(totva-meanc(totva'))^2;
nvoca=sqrt(sumc(nvocaa'));
dvocaa=totvsa^2;
dvoca=sqrt(sumc(dvocaa'));
ocva=100*(nvoca./dvoca);
ocva;
print " ";

/* check percentage of parameters within 95% conf int. */

print "check standard beta's percentage within 95% conf int";
cib=abs(bsum-meanc(bsum')).<(1.96*sesum);
sumc(cib')/(trial-1);
print " ";

print "check alternative beta's percentage within 95% conf int";
ciba=abs(bsuma-meanc(bsuma')).<(1.96*sesuma);
sumc(ciba')/(trial-1);
print " ";

print "check ols percentage within 95% conf int";
cio=abs(blsum-meanc(blsum')).<(1.96*selsum);
sumc(cio')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cit=abs(blsumo-meanc(blsumo')).<(1.96*selsumo);
sumc(cit')/(trial-1);
print " ";

print "check het ols percentage within 95% conf int";
cih=abs(bsumh-meanc(bsumh')).<(1.96*sesumh);
sumc(cih')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cil=abs(bsuml-meanc(bsuml')).<(1.96*sesuml);
sumc(cil')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cill=abs(blsumol-meanc(blsumol')).<(1.96*selsumol);
sumc(cill')/(trial-1);
print " ";

/* get information about percentage of significant parameters */
/* mean values */
print "check percentage of sig parms (standard beta)";
tstatbm=(abs(totm./totms)).>1.96;
tstatbv=(abs(totv./totvs)).>1.96;
sumc(tstatbm')/(trial-1);
sumc(tstatbv')/(trial-1);
print " ";

print "check percentage of sig parms (alternate beta)";
tstatba=(abs(bsuma./sesuma)).>1.96;
sumc(tstatba')/(trial-1);
print " ";

print "check percentage of sig parms (basic ols)";
tstatbo=(abs(blsum./selsum)).>1.96;
sumc(tstatbo')/(trial-1);
print " ";

print "check percentage of sig parms (logit trans ols)";
tstatbl=(abs(blsumo./selsumo)).>1.96;
sumc(tstatbl')/(trial-1);
print " ";

print "check percentage of sig parms (het ols)";
tstatboh=(abs(bsumh./sesumh)).>1.96;
sumc(tstatboh')/(trial-1);
print " ";

print "check percentage of sig parms (het ols - logit trans)";
tstatblh=(abs(bsuml./sesuml)).>1.96;
sumc(tstatblh')/(trial-1);
print " ";

print "check percentage of sig parms (ols -- log tranform)";
tstatbm=(abs(blsumol./selsumol)).>1.96;
sumc(tstatbm')/(trial-1);
print " ";


/* obtain information about y's */
print "descriptive stats for y's for successful runs";
meanc(ystats);
print "descriptive stats for y's for failed runs";
meanc(ystatf);

/* obtain information about x's */
print "descriptive stats for x's for successful runs";
meanc(xstats');
print "descriptive stats for x's for failed runs";
meanc(xstatf');

print "number of failed runs";
crash;
output off;